TK!2
=v
=u
=l
#1
:n
coeff
:v
.32432
.00565
.27196
.37294
.60025
.17577
.58748
.64068
.91559
.81142
.09197
E
#2
:n
Re
:v
E
#3
:n
Im
:v
E
=f
#1
:n
Bairstow
:c
Bairstow's algorithm for complex roots of real polynomial
:t
Procedure
:i
a,re,im
#1
:n
Bairstow
:z
" Notation:
" a    - name of list with polynomial cofficients in descending order
" re   - name of list with real components of the computed roots
" im   - name of list with imaginary components of the computed roots
" aux1,aux2  - names of scratch lists

" Description:
" Bairstow method for finding real and complex roots of real polynomials
" with repetitive extraction of quadratic factors using an internal
" Newton's iterative procedure.

if a[1]=0 then error:= highest_order_coefficient_cannot_be_zero
call blank(re)
call blank(im)
for n=length(a)-1 to 1 step -1
   if a[n+1]<>0 goto 1
   re[n]:= 0
   im[n]:= 0
next n
return
1: if n>1 goto 2
   (re[1],im[1]):= (-a[2]/a[1],0)     "Linear equation
   return
2: call listcopy(a,'aux1)
   n1:= 2*int((n+1)/2)
   'aux1[n+2]:= 0
   call blank('aux2)
   e1:= 1e-4
   index:= 1
A: if n1>2 goto B
     p:= 'aux1[2]/'aux1[1]     "Last quadratic equation
     q:= 'aux1[3]/'aux1[1]
     goto ST
B: p:= 1
   q:= 1
   for k=1 to 50                      "Iteration loop
C:    call listcopy('aux1,'aux2)
      for j=n1-2 to n1-4 step -2
         for i=1 to j+1
            'aux2[i+1]:= 'aux2[i+1] - p*'aux2[i]
            'aux2[i+2]:= 'aux2[i+2] - q*'aux2[i]
         next i
      next j
      r0:= 'aux2[n1+1]
      r1:= 'aux2[n1]
      s0:= 'aux2[n1-1]
      s1:= 'aux2[n1-2]
      v0:= -q*s1
      v1:= s0 - s1*p
      d0:= v1*s0 - v0*s1
      if abs(d0)>=1e-20 goto PQ
        p:= p + 5
        q:= q + 5
        goto C       " d0 too close to 0, recalculate for new p,q
PQ:  p1:= (s0*r1-s1*r0)/d0
     q1:= (r0*v1-v0*r1)/d0
     p:= p + p1
     q:= q + q1
     if min(max(abs(r0),abs(r1)),max(abs(p1),abs(q1)))<e1 goto ST
     if p*q=0 goto Loopend
     if max(abs(p1/p),abs(q1/q))<e1 goto ST
Loopend: next k
          error:= convergence_failure
ST:  s:= -p/2
     t:= s^2 - q
     if t<=0 goto Complex
       t:= sqrt(t)
     if or(mod(n,2)=0,abs(abs(s)/t-1)>e1) goto Real
       (re[index],im[index]):= (s+t*sgn(s),0)   "Single real root
       index:= index - 1
       goto Reset
Real:    (re[index],im[index]):= (s+t,0)        "Pair of real roots
         (re[index+1],im[index+1]):= (s-t,0)
         goto Reset
Complex: t:= sqrt(-t)                           "Pair of complex roots
         (re[index],im[index]):= (s,t)
         (re[index+1],im[index+1]):= (s,-t)
Reset: index:= index + 2
       for j=1 to n1-1
          'aux1[j+1]:= 'aux1[j+1] - p*'aux1[j]
          'aux1[j+2]:= 'aux1[j+2] - q*'aux1[j]
       next j
       n1:= n1-2
       if n1>0 goto A
       call delete('aux1)
       call delete('aux2)
/E
E
=r
#1
:r
"      * * *  Real and complex roots of a polynomial  * * *
:s
C
#3
:r
call Bairstow('coeff,'Re,'Im)
#5
:r
" This model exemplifies the usage of the procedure function 'Bairstow'
:s
C
#6
:r
" for finding all real and complex roots of a real-coefficient polynomial.
:s
C
#7
:r
" Type  = F >  to look at it.
:s
C
#9
:r
" An interactive table ROPE may be used for setting up the coefficients and
:s
C
#10
:r
" watching the root generation (type  = T > > F9 ).  The coefficients defined
:s
C
#11
:r
" in the list   coeff   represent a 10th order polynomial with two real and
:s
C
#12
:r
" eight complex roots (the solution takes a while if you have a slow machine
:s
C
#13
:r
" without a math co-processor).
:s
C
#15
:r
" Once the solution is found, press  F7  to display the roots in the complex
:s
C
#16
:r
" plane.
:s
C
=p
#1
:n
ROPE
:t
Line chart
:d
EGA
:i
"Roots of a Polynomial Equation in the Complex Plane "
:s
Yes
:z
Both
:g
No
:l
None
:h
"R e a l"
:v
"Imaginary"
:x
Re
:u
#1
:y
Im
:f
1
:p
Symbols
:c
"û"
E
=t
#1
:n
ROPE
:s
Screen
:v
Vertical
:f
1
:r
""
:c
""
:t
"Roots of a Polynomial Equation"
:u
#1
:l
coeff
:w
10
:h
""
#2
:l
Re
:w
10
:h
""
#3
:l
Im
:w
10
:h
""
E
=n
#1
:n

:d
15
:f

:^
Either Scientific or Decimal
:e
"."
:m
""
:p
None
:j
Left
:s
- Only
:q
0
:x
0
:l
""
:t
""
:v
"0"
=g
#1
:c
.000001
:i
10
:v
Yes
:l
60
:f
66
:m
0
:w
80
:b
Yes
:a
Yes
:r
Yes
:h
No
:t
1
:n
Yes
:d
LPT1
%Tr,f,0,1,0,1,0
